Graphic Solution of Differential Equation of the Second Order

SystemDiffeq[x1_,y1_,z1_,f_,g_,h_,n_]:=Module[{w=Table[0,{i,1,5},{j,1,n}]},w[[1,1]]=x1;w[[2,1]]=y1;w[[3,1]]=z1;Do[w[[4,i]]=f[w[[1,i]],w[[2,i]],w[[3,i]]];w[[5,i]]=g[w[[1,i]],w[[2,i]],w[[3,i]]];w[[1,i+1]]=w[[1,i]]+h;w[[2,i+1]]=w[[2,i]]+w[[4,i]]h;w[[3,i+1]]=w[[3,i]]+w[[5,i]]h,{i,1,n-1}]
;w];
f3[x_,y_,z_]:=z/(x+1);
g3[x_,y_,z_]:=((y-z)^2+x z)/(x+.1)^2;
f1[x_,y_,z_]:=z;
g1[x_,y_,z_]:=-x^2y-z;
f2[x_,y_,z_]:=z;
g2[x_,y_,z_]:=-0.1x +0.3y;
f4[x_,y_,z_]:=z;
g4[x_,y_,z_]:=-y+Exp[x];
f5[x_,y_,z_]:=z;
g5[x_,y_,z_]:=-1;
f6[x_,y_,z_]:=z;
g6[x_,y_,z_]:=-x;
f7[x_,y_,z_]:=z;
g7[x_,y_,z_]:=2x;
f8[x_,y_,z_]:=z;
g8[x_,y_,z_]:=1+y z;

pairs={{f1,g1},{f2,g2},{f3,g3},{f4,g4},{f5,g5},{f6,g6},{f7,g7},{f8,g8}};

Manipulate[With[{f=pairs[[pf,1]],g=pairs[[pf,2]]},With[{w=SystemDiffeq[x1,y1,z1,f,g,h,n]},With[{gr=Grid[{Join[{" "},Range[n]],Join[{Style["x",Italic]},w[[1,All]]],Join[{Style["y",Italic]},w[[2,All]]],Join[{Style["z",Italic]},w[[3,All]]],Join[{Style["y",Italic]'},w[[4,All]]],Join[{Style["z",Italic]'},w[[5,All]]]},Frame->All]},If[shg,Text@Column[{TraditionalForm[y'==f[x,y,z]],TraditionalForm[z'==g[x,y,z]],Show[If[sn,With[{s=Quiet[NDSolve[{y''[x]==g[x,y[x],y'[x]],y[x1]==y1,y'[x1]==z1},y,{x,x1,x1+(n-1)h}]]},Plot[Evaluate[{y[x],y'[x]}/.s],{x,x1,x1+h(n-1)},PlotStyle->Automatic,PlotRange->{{-.1,1.3},{-0.2,1.6}}]],Graphics[{}]],Graphics[{PointSize[0.02],{RGBColor[0, 0.866667, 0.866667],Thickness[0.008],Line[Table[{w[[1,i]],w[[2,i]]},{i,1,n}]]},{Green,Thickness[0.008],Line[Table[{w[[1,i]],w[[3,i]]},{i,1,n}]]},If[shexp,{Red,Line[{{w[[1,n-1]],w[[2,n-1]]},{w[[1,n]],w[[2,n-1]]}}],Arrowheads[Small],Arrow[{{w[[1,n]],w[[2,n-1]]},{w[[1,n]],w[[2,n]]}}],Line[{{w[[1,n-1]],w[[3,n-1]]},{w[[1,n]],w[[3,n-1]]}}],Arrow[{{w[[1,n]],w[[3,n-1]]},{w[[1,n]],w[[3,n]]}}],Text[Style["h",Italic,10],{(w[[1,n-1]]+w[[1,n]])/2,w[[2,n-1]]},{0,1}],Text[Style[Row[{Style["h",Italic]," ",TraditionalForm[y'[(Subscript[x, n-1])]]}],10],{w[[1,n]],(w[[2,n-1]]+w[[2,n]])/2},2{-1,0}],
Text[Style[Row[{TraditionalForm[(Subscript[y, n-1])]}],10],{w[[1,n]],w[[2,n-1]]},3{-1,0}],
Text[Style[Row[{TraditionalForm[(Subscript[y, n])]}],10],{w[[1,n]],w[[2,n]]},3{-1,0}],Text[Style[Row[{Style["h",Italic]," ",TraditionalForm[z'[(Subscript[x, n-1])]]}],10],{w[[1,n]],(w[[3,n-1]]+w[[3,n]])/2},2{-1,0}],Text[Style[Row[{TraditionalForm[(Subscript[z, n-1])]}],10],{w[[1,n]],w[[3,n-1]]},3{-1,0}],Text[Style[Row[{TraditionalForm[(Subscript[z, n])]}],10],{w[[1,n]],w[[3,n]]},3{-1,0}]},{}],Table[Text[Style[(Subscript[Style["x",Italic], i]),10],{w[[1,i]],0},3{0,1}],{i,1,n}],Point[Table[{w[[1,i]],w[[2,i]]},{i,1,n}]],Point[Table[{w[[1,i]],w[[3,i]]},{i,1,n}]],Table[Tooltip[Point[{w[[1,i]],w[[2,i]]}],Style[Row[{(Subscript["P", i]),"(",(Subscript[Style["x",Italic], i]),", ",(Subscript[Style["y",Italic], i]),")"}],10]],{i,1,n}],Table[Tooltip[Point[{w[[1,i]],w[[3,i]]}],Style[Row[{(Subscript["T", i]),"(",Style[(Subscript[Style["x",Italic], i]),10],", ",Style[(Subscript[z, i]),10],")"}],10]],{i,1,n}]}],AspectRatio->Automatic,ImageSize->{400,400},PlotRange->{{-.1,1.3},{-0.2,1.6}},Axes->True,GridLines->{Table[x1+i h,{i,0,n-1}],None},GridLinesStyle->GrayLevel[0.8],AxesOrigin->{0,0}]}],Text@Column[{" ",TraditionalForm[y'==f[x,y,z]],TraditionalForm[z'==g[x,y,z]]," ",gr}]]]]],

{{pf,1,Row[{"(",Style["f",Italic],",",Style["g",Italic],")"}]},Range[8]}," ",
"initial values",
{{x1,0.1,TraditionalForm[(Subscript[x, 1])]},0.1,0.2,0.01,Appearance->"Labeled", ImageSize-> Tiny},{{y1,.3,TraditionalForm[(Subscript[y, 1])]},.3,.6,0.01,Appearance->"Labeled", ImageSize-> Tiny},{{z1,1.4,TraditionalForm[(Subscript[z, 1])]},0.3,1.5,0.01,Appearance->"Labeled", ImageSize-> Tiny}," ",
"step",
{{h,0.2,Style["h",Italic]},0.1,0.2,0.001,Appearance->"Labeled", ImageSize-> Tiny}," ","number of points",
{{n,3,Style["n",Italic]},2,6,1,Appearance->"Labeled", ImageSize-> Tiny}," ","graphics/calculations",
{{shg,True,""},{False,True}}," ","show 'exact' solution",{{sn,False,""},{False,True}}," ","show explanation",{{shexp,False,""},{False,True}},
ContentSize->{430,490},
SaveDefinitions->True]

This Demonstration shows the Euler–Cauchy method for the graphics solution of the initial problem of a second-order differential equation. An example of such an equation is y''+y'+(x^2)y=0, with derivatives from now on always taken with respect to x. This equation can be written as a pair of first-order equations, y'=z, z'=-(x^2)y+z.

More generally, the method to be described works for any system of two first-order differential equations y'=f(x,y,z), z'=g(x,y,z) with initial conditions y((Subscript[x, 1]))=(Subscript[y, 1]), z((Subscript[x, 1]))=(Subscript[z, 1]). The particluar kind of systems used as examples here, y''=g(x,y,y'), are of that general type if we introduce z=y' to get the system y'=z, z'=g(x,y,z).

The method consists of simultaneously calculating approximations of y(x) (cyan) and z(x) (green):

(Subscript[y, i+1])=(Subscript[y, i])+h y'((Subscript[x, i]))=(Subscript[y, i])+h f((Subscript[x, i]),(Subscript[y, i]),(Subscript[z, i])),

(Subscript[z, i+1])=(Subscript[z, i])+h z'((Subscript[x, i]))=(Subscript[z, i])+h g((Subscript[x, i]),(Subscript[y, i]),(Subscript[z, i])), i=1, ..., n-1.

The pairs ((Subscript[x, i]),(Subscript[y, i])) are the coordinates of points (Subscript[P, 1]), (Subscript[P, 2]), ..., (Subscript[P, n]) that form the so-called Euler's polygonal line that approximates the graph of the function y(x). In the same way, the pairs ((Subscript[x, i]),(Subscript[z, i])) are the coordinates of points (Subscript[T, 1]), (Subscript[T, 2]), ..., (Subscript[T, n]) that form Euler's polygonal line that approximates the graph of the function z(x).

The Euler method is the most basic approximation method. We compare it with more advanced methods given by the built-in Mathematica's function NDSolve.

V. I. Smirnoff, Lectures in Higher Mathematics (in Russian), Vol. II, Nauka, Moscow 1967, p. 50.

L. Euler, Institutionum calculi integralis volumen primum, section 2, chapter 7: De integratione aequationum differentialium per approximationem, written 1763, published 1768.

Contributed by: Izidor Hafner